Spatio-temporal dynamics of bumblebees foraging under predation risk 
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We analyze 3D flight paths of bumblebees searching for nectar in a laboratory experiment with and without 
predation risk from artificial spiders. For the flight velocities we find mixed probability distributions reflecting 
the access to the food sources while the threat posed by the spiders shows up only in the velocity correlations. 
The bumblebees thus adjust their flight patterns spatially to the environment and temporally to predation risk. 
Key information on response to environmental changes is contained in temporal correlation functions, as we 
explain by a simple emergent model. 
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Quantifying foraging behavior of organisms by statistical 
analysis has raised the question whether biologically relevant 
search strategies can be identified by mathematical model- 
ing OJ-0]. For sparsely, randomly distributed, replenishing 
food sources, the Levy flight hypothesis predicts that a ran- 
dom search with jump lengths following a power law mini- 
mizes the search time 0-ISD- Experimental evidence 1 13-13] 
and further theoretical analyses [l5ll supporting this hy- 
pothesis were challenged by refined statistical data analyses 
[16 -IH and more detailed theoretical modeling j6],[2(J[2l|]. A 
crucial problem is how dispositions of a forager like memory 
1 2211 or senso ry p erception 112311 . as well as properties of the 



environment IU2 U L3 U 24 1 27I1. can be tested in a statistical for- 
aging analysis llRji |5l |7|] . Especially for data obtained from 
foraging experiments in the wild, it is typically not clear to 
which extent extracted search patterns are determined by for- 
ager dispositions, or reflect an adjustment of the dynamics of 
organisms to the distribution of food sources and the presence 
of predators [T3I1 . This problem can be addressed by 

statistically quantifying search behavior in laboratory exper- 
iments where f oraging conditions are varied in a fully con- 




trolled manner II 1 3L 12411 . Such an experiment has been per- 



formed by Ings and Chittka [28, 29], who studied the foraging 
behavior of bumblebees with and without different types of 
artificial spiders mimicking predators. 

Here we ask the question whether changes of environmental 
conditions as performed in the experiment by Ings and Chittka 
lead to changes in the foraging process. We answer this ques- 
tion by a statistical analysis of the bumblebee flights recorded 
in this experiment on both spatial and temporal scales. For 
this purpose, we extract both flight velocity probability distri- 
butions and temporal velocity autocorrelation functions from 
the data. Surprisingly, we find that the crucial quantity to un- 
derstand changes in the bumblebee dynamics under predation 
risk is not the velocity distribution but the velocity correlation 
function, which reveals non-trivial dynamics on different time 
scales. We reproduce these changes by a simple Langevin 
equation modeling a repulsive interaction between insect and 
predator. In order to construct mathematical models repro- 
ducing the foraging of organisms that interact with the envi- 
ronment, our results suggest to shift the focus from scale-free 



Figure 1. Diagram of the foraging arena together with part of the 
flight trajectory of a single bumblebee. The bumblebees forage on 
a grid of artificial flowers on one wall of the box. While being on 
the landing platforms, the bumblebees have access to food supply. 
All flowers can be equipped with spider models and trapping mech- 
anisms simulating predation attempts. 



approaches QUI OH to the statistical quantification of spatio- 
temporal changes in the foraging dynamics. 



In the experiment [28] bumblebees (Bombus terrestris) 
were flying in a cubic arena of w 75 cm side length by for- 
aging on a 4x4 vertical grid of artificial yellow flowers on 
one wall. The 3D flight trajectories of 30 bumblebees, tested 
sequentially and individually, were tracked by two high frame 
rate cameras (At = 0.02 s). On the landing platform of each 
flower, nectar was given to the bumblebees and replenished 
after consumption. The short trajectory in Fig.[T]shows a typ- 
ical flight path of a bumblebee foraging in the arena. To an- 
alyze differences in the foraging behavior of the bumblebees 
under threat of predation, artificial spiders were introduced. 
The experiment was staged into three phases: (1) spider-free 
foraging, (2) foraging under predation risk and (3) a mem- 
ory test one day later. Before and directly after stage (2) the 
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zero mean: p a {v) = N a (v), (d) Mixture of two normal dis- 
tributions: /5 a ,(Ti.cr 2 («) = aN ai (v) + (1 — a)Ncr 2 (v), where 
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Figure 2. Estimated velocity distributions (main part) and Quantile- 
Quantile probability plot of a Gaussian mixture as the best fit (in- 
set). Semi-logarithmic plot of the normalized histogram of velocities 
v y parallel to the y-axis in Fig. [T] (black crosses) for a single bum- 
blebee in the spider-free stage (1) together with a Gaussian mixture 
(red line), exponential (blue dotted), power law (green dashed), and 
Gaussian distribution (violet dotted), fitted via maximum likelihood 
estimation. The inset shows quantiles of v y (in m /s) of a single bum- 
blebee against quantiles of an estimated mixture of two Gaussians. 
An ideal match would yield a straight line. The dashed red lines 
show 20 surrogate data sets of the same size. 



bumblebees were trained to forage in the presence of artificial 
spiders, which were randomly placed on 25% of the flowers. 
A spider was emulated by a spider model on the flower and 
a trapping mechanism which held the bumblebee for two sec- 
onds to simulate a predation attempt. In (2) and (3) the spiders 
models were present but the traps were inactive in order to an- 
alyze the influence of previous experience with predation risk 
on the bumblebees' flight dynamics. To determine whether 
the detectability of the spiders is an important factor, half of 
the bumblebees were trained on easily visible (white) spider 
models and half of them on yellow models, which meant that 
spiders were camouflaged on the yellow flowers; see Ings and 
Chittka i28ll for further details of the experiment. 

Figure |2] shows a typical normalized histogram of the hori- 
zontal velocities parallel to the flower wall (cf. y-direction in 
Fig. [TJ for a single bumblebee. The histograms are charac- 
terized by a peak at low velocities and vary in the different 
spatial directions due to asymmetries induced by physical and 
biological constraints as well as the spatial arrangement of the 
flowers. Direct fitting of distributions on the histogram and a 
visual comparison with some assumed distribution was shown 
to be unreliable II 1711 . as is illustrated by Fig. [2] only the power 
law and the Gaussian distribution can be ruled out by visual 
inspection. However, the Gaussian mixture and an exponen- 
tial function appear to be equally likely. Therefore we use 
the maximum likelihood method for a number of candidate 
distributions to obtain the optimal parameters for each candi- 
date and then compare the different distribution types by their 
weights using the Akaike information criterion 1 1611 . Our can- 
didate distributions are: (a) Exponential: p\(v) = ce~ x > v \, 



= 1, 2, and < a < 1. Details of 



(b) Power law: p^(v) 
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, (c) Normal distribution with 



this analysis are described in the Supplemental Material 04 1. 

For the data sets of all bumblebees and in all stages of the 
experiment the Akaike weights show that a mixture of two 
Gaussians is the preferred distribution of the tested candidates 
(see Table I in the Supplemental Material 13411 ). However, 
they do not inform us if the best of the candidates is actually 
a good model: if all of the candidates are far off the real dis- 
tribution, the Akaike weights could highlight one of them as 
the best of the poor fits. As a supplementary qualitative test to 
which extent the estimated distribution with the largest Akaike 
weight deviates from the data over the whole range variables, 
we use Quantile-Quantile (Q-Q) probability plots. The inset 
of Fig.|2]shows the Q-Q plot of the mixture of two Gaussians 
against the experimental data of a single bumblebee and 20 
surrogate data sets. Each of the surrogate data sets consists of 
independently identically distributed random numbers drawn 
from the estimated Gaussian mixture and has the same number 
of data points as the real data for comparison with statistical 
fluctuations. The Q-Q plot shows that the deviations of the 
experimental data from the mixture of two Gaussians is not 
larger than the expected deviations due to the finite quantity 
of data. 

The Gaussian mixture for the velocities is generated by 
different flight behavior near a flower versus in open space, 
which bears some resemblance to intermittent dynamics lalTi 
2511 . This has been verified by splitting the data into flights 
far from the flower wall vs. flights in the feeding zone. The 
latter was defined by a cube of side length 9 cm around each 
flower in which the velocities are determined by approach- 
ing a flower and hovering behavior. This separation of differ- 
ent flight phases is thus adapted to accessing the food sources 
and explains the origin of Gaussian distributions with different 
variances in both spatial regions. Because of the absence of a 
sparse distribution of food sources, there was no reason to ex- 
pect Levy-type probability distributions Jgl. Surprisingly, by 
comparing the best fits to these distributions for the different 
stages of the experiment, we could not detect any differences 
in the velocity distributions between the spider-free stage and 
the stages where artificial spider models were present, as is 
shown in Table II of the Supplemental Material 1 34] . The 
parameters of the Gaussian mixture vary between individual 
bumblebees but there is no systematic change due to the pres- 
ence of predators. 

Hence, we examined the velocity autocorrelations for com- 
plete flights from flower to flower. The autocorrelations 
have been computed by averaging over all bumblebees while 
weighting with the amount of data available for each time in- 
terval; see the Supplemental Material for details 13411 . Figure[3] 
shows the velocity autocorrelations in the x- and y-directions 
for different stages of the experiment. In the x-direction per- 
pendicular to the wall the velocities are always anti-correlated 
for times around 0.5 s (Fig.Oa)), which is due to the tendency 
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of the bumblebees to quickly return to the flower wall. How- 
ever, the flights with long durations between flower visits be- 
come more frequent for stages (2) and (3) where the bumble- 
bees were exposed to predation risk compared with stage (1) 
(inset of Fig.|3ja)). This is also reflected in a small shift of the 
global minimum in the correlations for stages (2) and (3) away 
from the origin. The v x -autocorrelations thus display similar 
functional forms but with quantitative changes between the 
different stages. In contrast, for the v y -autocorrelations the 
functional forms change profoundly: Parallel to the flower 
wall the velocities are anti-correlated in the presence of spi- 
ders for 0.7 s < t < 2.8 s, while for the spider-free stage the 
correlations remain positive up to 1.7 s (Fig. Ob)). The verti- 
cal z-direction is similar to the y-direction with a weaker de- 
pendence on the presence of predators. As the limited amount 
of data causes variations in the autocorrelations between in- 
dividual bumblebees, we resampled the result by leaving the 
data of each single bumblebee out (jackknifing). The resam- 
pling (inset of Fig. Ob)) confirms that the positive autocorre- 
lations of v y are not a numerical artifact. 

Our statistical analysis of the experimental data has thus 
revealed that differences in the foraging behavior of bumble- 
bees, triggered by predation risk, show up in changes of the 
velocity autocorrelation functions only, and not in modifica- 
tions of the velocity probability distributions. These changes 
are consistent with a more careful search: When no threat of 
predators is present, the bumblebees forage more systemati- 
cally with more or less direct flights from flower to flower, 
arching away from the flower wall. Under threat the trajecto- 
ries become longer and the bumblebees change their direction 
more often in their search for food sources, rejecting flowers 
with spiders, as is supported by Fig. [3] Further analysis rules 
out that the main features of the correlation functions are in- 
duced by the geometry of the experiment: In Fig. 0a), all 
flight time distributions display maxima around Tj « 0.5 s 
suggesting that times below ~ 2 s are primarily related to 
flights between flowers. Boundary effects are only evident 
for flight times that fall within the tail of the distributions. 
The anti-correlations in the y- and z-directions parallel to the 
flower wall thus cannot be induced by the walls but are gen- 
erated by a reversal of directions at flowers under predatory 
threat. For the x-direction, the return to the flower wall is re- 
sponsible for the anti-correlation at small delay times, not the 
opposite wall, which is too far away to have a significant ef- 
fect. Another interesting aspect is that Ings and Chittka 112 811 
report that bumblebees increase the time they spend inspect- 
ing flowers bearing camouflaged spiders compared to conspic- 
uous ones. However, we did not detect any change in veloc- 
ity distributions or autocorrelations in this case, which sug- 
gests that the bumblebees perform longer localized inspection 
flights without changing their velocities. 

In order to understand the changes shown in Fig. [3] we 
model the dynamics of v y by the Langevin equation 
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Figure 3. Autocorrelation of the velocities at different experimen- 
tal stages: without spiders (red triangles), under threat of predation 
(green circles), and under threat a day after the last encounter with 
the spiders (blue crosses), (a) In the x-direction perpendicular to the 
wall the velocities are anti-correlated for small times (~ 0.5 s) due 
to short flights from one flower to a nearby flower back at the flower 
wall. Inset: the distribution of flight times Tj for each stage shows 
a corresponding maximum for these short jumps. Under threat of 
predation (dotted) long flights become more frequent, (b) The cor- 
relation of v y parallel to the wall shows the effect of the presence of 
spiders on the flight behavior of the bumblebees. The inset shows 
the resampled autocorrelation for the spider-free stage in the region 
where the correlation differs from the stages with spider models, 
which confirms that the positive autocorrelations are not a numeri- 
cal artifact. 



where 77 is a friction coefficient and £ Gaussian white noise. 
The potential U mimics an interaction between bumblebee 
and spider. Specific data analysis shows that this force is re- 
pulsive and dominates any hovering behavior in the velocity 
correlation decay. Computer simulations of the above equa- 
tion reproduce a change from positive to anti-correlations by 
increasing the repulsive force. Details are discussed in the 
Supplemental Material l34tl . Note that the correlation decay 
displayed in Fig.[3]rules out a mathematical modeling in terms 
of ordinary correlated random walks or Levy walks, which 
predict velocity correlations to decay strictly exponentially 
12, 26] or algebraically |3(| 31], respectively. 



We emphasize that the experiment analyzed in this Letter 
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does not match the conditions of the Levy flight hypothesis 
J9t]. Levy flights and Levy walks predict scale-free probabil- 
ity distributions [7] and gene rate trivial functional forms for 
the velocity correlations |3(| 31]. Accordingly, experiments 
testing this hypothesis have focused on probability distribu- 
tions, not on correlation decay [13-13]. However, our results 
demonstrate that velocity autocorrelations can contain crucial 
information for understanding foraging dynamics, here in the 
form of a highly non-trivial correlation decay emerging from 
an interaction between forager and predator. Identifying such 
an emergent property in contrast to adaptive behavior, as we 
do with our simple model, has been highlighted as a crucial 
problem in foraging dynamics 111 311 . In addition, we observe 
a spatial variation of the velocity distributions. These find- 
ings illustrate the presence of different flight modes governing 
the foraging dynamics on different scales of time and space. 
Our results thus indicate that taking scale-free distributions 
as a paradigm beyond the conditions of validity of the Levy 
flight hypothesis might be too restrictive an approach in order 
to capture complex foraging dynamics. A variety of mech- 
anisms may naturally lead to different foraging dynamics on 
different length and time scales, e.g., individuality of animals 
I IH 32, 33], an intermittent switching between quasi-ballistic 
persistent dynamics and localized search modes |0, flil]. o r 
quantities over which one has averaged like time of day 1 1311 . 
As ignoring these mechanisms can lead to spurious power 
laws 1 3, 171. it is important to look for the reasons of the 
occurrence of non-trivial distributions like mixtures, e.g., an- 
imals switching between different search modes. These mix- 
tures may not always be optimal distributions for a particular 
search problem, but they are easy to produce, composable and 
flexible enough such that differences to some optimal distri- 
bution might not be large enough to give rise to evolutionary 
pressure |4J]. 

In summary, the fundamental question 'What is the math- 
ematically most efficient search strategy of foraging organ- 
isms?' has, under specific conditions 12111 . been answered by 
the Levy flight hypothesis Jil0]. This question is well-posed 
under precise foraging conditions and has the big advantage 
that it is amenable to mathematical analysis. However, it does 
not capture the full complexity of a biological foraging prob- 
lem QV|] , which incorporates both the dependence of foraging 
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on 'internal' conditions of a fora ger ( sensory perception 12311 . 



memory 12211 . individuality 0191 1321 |33[1 ) as well as 'exter 



nal' environmental constraints (distribution of food sources 
lHHlll^], day-night cycle Q, predators JUB). Asking 
about the range of applicability of the Levy flight hypothe- 
sis leads to the over-arching question 'How can we statisti- 
cally quantify changes in foraging dynamics due to interac- 
tions with the environment?', which requires to identify suit- 
able measurable quantities characterizing such changes. This 
question highlights the need to better understand, and more 
carefully analyze, the interplay between forager and environ- 
ment, which will yield crucial information for constructing 
more general mathematical foraging models. 

We thank Holger Kantz for valuable support and Nicholas 
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SUPPLEMENTAL MATERIAL 

1. Statistical data analysis 

To capture bumblebee flights only, we exclude any crawling 
behavior on the landing platforms by also removing all data 
within a 1 cm boundary region of each platform. The size of 
this boundary is based on the size of the bumblebees, which 
have a height of approximately 1 cm. While smaller cutoffs 
would not exclude all crawling behavior, the cutoff can be in- 
creased robustly within reasonable bounds. We have checked 
that, e.g. a 2 cm cutoff does not have any influence on any of 
the analyzed quantities, as the amount of the data which would 
be excluded in addition is very small. This leaves from 2000 
to 15000 data points (average: 6000) per bumblebee for each 
stage. We select the best model for the velocity distributions 
by maximum likelihood estimation and Akaike and Bayesian 
weights for our candidate distributions [16] for |w| > 2.5 cm /s. 
Given a set of measured velocities D = {i>x, v%, v n } and 
a probability density function p\(v), where A is a vector of k 
parameters, the log-likelihood of the probability density func- 
tion for a finite resolution of the data (Av = 5 cm /s) simplifies 
to 



/max(b) 
nin(b) 



p\(v)dv 



bins 



where h(b) is the observed frequency in bin b. 

For each candidate distribution p\.,i <G {1, 2, 3}, we max- 
imize the log-likelihood In Li w.r.t. A^ locally with a Nelder- 
Mead algorithm by using a Monte Carlo method to find the 
global maximum. To find the preference between the differ- 
ent model distributions whose likelihoods Li are maximized 



at Wl iax the information criteria are 

Id = -2\n{L t {\? ax \D)) + S (n)h 

with s(n) = 2 for the Akaike information criterion and 
s(n) = ln(n) for the Bayesian information criterion as a 
penalty on the number of parameters fcj. The best model, de- 
noted by *, is the one which minimizes the information cri- 
terion 7C* = mm(ICi). The Akaike/Bayesian weights then 

i 

give the preference of each model over the others as a proba- 
bility 



W: 



-(IC'i-IC,)/2 



where a normalizes the weights to ^\ Wi = 1. 

The choice of the information criterion makes no strong dif- 
ference for the model selection in this experiment. With the 
Akaike information criterion the Gaussian mixture is chosen 
with a weight of over 95% for all bumblebees and all experi- 
mental stages. The Bayesian information criterion agrees with 
the Akaike information criterion on 90% of all data sets. For 
the other 10% it prefers a single Gaussian or an exponential 
distribution - these data sets turned out to be those with the 
least amount of data available. 

To compute the autocorrelation function v ac {r) of the flight 
velocities 



v ac (r) 



((«(*)- fi)(v(t + t)-h)) 



we average over all bumblebees and over time in all flights 
that are complete from starting on one flower to landing on 
the next. We exclude flights containing gaps and correlation 
terms, where in-between time t and f + ra flower was visited. 



Table I. Model weights and estimated parameters. Akaike and Bayesian weights both give preference to the mixture of two Gaussians for 
v y for most of the bumblebees. The weights are estimated individually and their mean and standard deviation (in brackets) are shown. The 
distribution parameters are also estimated individually for each bumblebee in each stage. 

Model: (a) Exponential (b) Power law (c) Gaussian (d) Gaussian Mixture 



Akaike weight 
Bayesian weight 
Parameters 
average (bumblebees) 
stddev (bumblebees) 



0.00 (0.00) 
0.04(0.18) 
A 

5.61 

1.07 



0.00 (0.00) 
0.00 (0.00) 

/' 
1.11 
0.16 



0.04(0.19) 
0.08 (0.26) 
a 

0.25 

0.03 



</ 
0.67 
0.13 



0.96 (0.19) 
0.88 (0.30) 

0.06 
0.04 



0.29 
0.03 



Table II. Weights and estimated parameters of the Gaussian mixture for the different experimental stages. Weights and parameters are estimated 
for each bumblebee. Shown are the mean over all individuals and the standard deviation (in brackets). The mixture of two Gaussians is the 
best fit in all stages. In the parameters of the distribution we observe no significant effect of the threat of predators on the bumblebees. 



Stages 


Akaike weight 


Bayesian weight 


a 


ai 


CT2 


(1) Without spiders 


0.97 (0.15) 


0.93 (0.23) 


0.64(0.11) 


0.06 (0.02) 


0.29 (0.03) 


(2) Under predation risk 


0.99 (0.04) 


0.90 (0.27) 


0.68 (0.13) 


0.06 (0.02) 


0.29 (0.02) 


(3) With risk, 1 day later 


0.89 (0.29) 


0.80 (0.38) 


0.72 (0.16) 


0.07 (0.07) 


0.30 (0.03) 



2 



App{x reh y re |) 



0.03 r 

0.02 
0.01 - 


-0.01 
-0.02 
-0.03 




0.08 



0.06 



-0.06 



Figure SI. Predator avoidance of bumblebees at flowers, Eq. Q}, ex- 
tracted from the experimental data. Hovering behavior in front of a 
flower is represented by the positive spike directly at the flower cen- 
ter, while the negative region behind this spike reflects the avoidance 
in the flights towards a flower. 

2. Mathematical modeling of bumblebee foraging 

The effect of the presence of a spider on the probability 
of a bumblebee to fly in front of a flower can be measured 
by computing the difference between the position densities at 
stage (1) and (2) as a function of the positions parallel to and 
near {x <5 cm) the flower wall, 

Ap p (y re i,Z re i) = pp(Vrel,Zrel) ~ Py illrel , Z re i) , (1) 

where the positions (y re i , Zrei ) are relative to the nearest 
flower center. This predator avoidance extracted from the ex- 
perimental data is shown in Fig. [Si] Two different types of 
behavior can be seen: First, there is a small increase in the 
amount of hovering, i.e. inspection flights near the flower plat- 
form when a spider model is present 13,12], which is consis- 
tent with Ref. 1 3fo. However, more important is the local min- 
imum representing the avoidance of flowers infected by spi- 
ders. This effect is strongest 3 cm above the dangerous flow- 
ers, because the flowers are predominantly approached from 
above. The avoidance behavior affects not only flights near 
the flower wall but can still be detected further away from it. 
Comparing dangerous and safe flowers at stage (2) only con- 
firms that avoidance is the dominant effect for search flights. 

The avoidance of spider-infected flowers together with the 
spatial switching of flight modes discussed in the main part of 
our Letter can be modeled by the Langevin Equation 

: (t)=v(t) 



dt 
dv 



{t) = -nv(t)-VU(r(t))+Z(T,t), 



(2) 



where r\ is a friction coefficient and £ white Gaussian noise 
with standard deviation depending on the flight mode as 
a function of the position, £(r,i) = Xfz( r )£iM + (1 — 
Xfz( r ))&(t)- Here r = (x, y, z) T is the position of the bum- 
blebee at time t, Xf z (r) is the indicator function of the feeding 



zone, which is equal to one whenever the bumblebee is in the 
cube around a flower as defined before, and £j , i = 1, 2 is 
Gaussian noise with two different variances. The potential U 
models an interaction between bumblebee and spider in form 
of a repulsive force exerted by the spider onto the bumblebee, 
for which we assume that the potential maxima are located 
near infected flowers. 

When the mechanism generating the correlation functions 
shown in Fig. 3 is not the focus of the investigation, it suf- 
fices to consider a reduced version of Eqs. (|2]) in form of the 
effective Langevin equation 



= Xfz(r)Ci(*) + (l-Xfz(r))Ca(*). 



(3) 



This equation describes the spatially varying hovering and 
search modes by using noise Q , i = 1, 2., which models the 
impact of the potential U together with the noise £. Further 
data analysis shows that excluding hovering has no significant 
impact on the velocity autocorrelations, which are dominated 
by the search flights. This is in full agreement with Fig. 3, 
where the time scale for the predator-induced anti-correlation 
(Fig. 3(b)) is larger than the time scale for flights between 
neighbouring flowers (Fig. 3(a)). Hence, we model d(t) as a 
vector of Gaussian white noise with the smaller variance a\ 
given in Table U which describes the hovering. The search 
flights from flower to flower are reproduced by the correlated 
Gaussian noise vector with variance cr| and the autocor- 
relations vf c (r) , i — x,y shown in Fig. 3. The advantage of 
this model is that it is directly based on our data analysis. 

We now focus on the different aspect of understanding the 
biophysical mechanism that generates the anti-correlations of 
the velocities parallel to y shown in Fig. 3(b). Starting from 
the full model Eqs. (0, it suffices to select the search mode 
only by setting £(r, i) = ^(t) thus neglecting any spatial 
variations of the noise. This yields the Langevin equation 



dt 



(t) = -r}v y {t) 



dU 
dy 



(4) 



which was already stated in the main part as the main equa- 
tion. A rough approximation for the repulsive force is pro- 
vided by a periodic potential with maxima at dangerous flow- 
ers, 



U(v) 



It COS 27T 



y_ 
yo 



(5) 



where yo is the mean distance between spiders and u the 
strength of the repulsion. 

We integrated this Langevin equation via an Euler- 
Maruyama method under variation of u by computing the au- 
tocorrelation function Vy C of the generated data. Figure [S2] 
shows Vy C by increasing the repulsion strength u. The cor- 
relation function changes from positive correlations to anti- 
correlations in a range of delay times t comparable to the 
changes in the correlation function of the experimental data of 
Fig. 3(b). This qualitatively reproduces our experimental find- 
ings from first principles. Note that the oscillations for higher 
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Figure S2. Autocorrelation function of the velocities v y for the 
Langevin model Eqs. 10, lO modeling predation threat by different 
strengths of a repulsive potential. Shown are results from computer 
simulations without (u — 0; red triangles, upper line) and with pre- 
dation threat (u — 0.5m 2 /s 2 ; green circles, lower line). These re- 
sults should be qualitatively compared with the experimental findings 
Fig. 3(b). 



t in Fig. [S2] would be suppressed in a higher-dimensional 
model. The other directions can be treated analogously, e.g., 
by including an x-dependent term in the potential for the at- 
traction of the bumblebees to the flower wall. A stochastic 



analysis of Langevin equations with periodic potentials can 
be found, e.g., in Ref. |4J]. The effect of the harmonic poten- 
tial on the creation of negative velocity correlations can also 
be calculated analytically JsJ] . 

We emphasize that our model Eqs. (HI,© provides only a 
qualitative description of the biophysical mechanism generat- 
ing the change in the correlations of the bumblebee velocities 
under predation threat. For a quantitative comparison to the 
experimental data a much more detailed model would be nec- 
essary, which needs to include the random positioning of the 
spiders and the general attractive force exerted by the flowers 
onto the bumblebees. Modeling the three-dimensional nature 
of the potential would also be important: Notice, e.g., the lo- 
cal maximum of Vy C around r ~ 2.5 which is an artifact of 
the one-dimensional modeling of spider avoidance. However, 
as it is difficult to reliably estimate the parameters of the po- 
tential, such a quantitative comparison is beyond the scope of 
our Letter. 
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